#This to obtain:
#table_A3.tex

rm(list = ls())
library("broom")
library("dplyr")
library("grid")
library("tidyr")
library("sandwich")
library("lfe")
library("broom")
library("dplyr")
library("tidyr")
library("grid")
library("stargazer")
library("lmtest")
library("multiwayvcov")
library("ggplot2")
library("glue")
library("spdep")
library("stringr")
library("lemon")
library("readxl")
library("spatialreg")

#setwd("C:/Users/bogdanp/")
setwd("/Users/bgpopescu/")

#Reading Polygons
HRV_adm0 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm0_wgs")

HRV_adm3 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm3")
HRV_adm3<-st_simplify(HRV_adm3, dTolerance = 0.005)
HRV_adm3<-subset(HRV_adm3, select = c(ID_2))
HRV_adm3_data = read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)

final_sp<-left_join(HRV_adm3, HRV_adm3_data, by = c("ID_2"="ID_2"))
data<-subset(final_sp, bosnia_border_mun==0)


########
#Models#
########
lm_elevation<-lm(elevation~treat+
                         POINT_X  + POINT_Y + 
                         zagreb_distance +  
                         bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                         quad1 + quad2 + quad3, data=data)

lm_slope<-lm(slope~treat+
                   POINT_X  + POINT_Y + 
                   zagreb_distance +  
                   bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                   quad1 + quad2 + quad3, data=data)

lm_tmp_mean_2010<-lm(tmp_mean_2010~treat+
               POINT_X  + POINT_Y + 
               zagreb_distance +  
               bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
               quad1 + quad2 + quad3, data=data)

lm_pre_mean_2010<-lm(pre_mean_2010~treat+
                       POINT_X  + POINT_Y + 
                       zagreb_distance +  
                       bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                       quad1 + quad2 + quad3, data=data)

lm_maizemed<-lm(maizemed~treat+
                       POINT_X  + POINT_Y + 
                       zagreb_distance +  
                       bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                       quad1 + quad2 + quad3, data=data)

lm_river_total<-lm(river_total~treat+
                  POINT_X  + POINT_Y + 
                  zagreb_distance +  
                  bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                  quad1 + quad2 + quad3, data=data)

lm_river_density<-lm(river_density~treat+
                     POINT_X  + POINT_Y + 
                     zagreb_distance +  
                     bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                     quad1 + quad2 + quad3, data=data)

lm_trade_route_1450_density<-lm(trade_route_1450_density~treat+
                       POINT_X  + POINT_Y + 
                       zagreb_distance +  
                       bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                       quad1 + quad2 + quad3, data=data)


lm_trade_center_count<-lm(trade_center_count~treat+
                                  POINT_X  + POINT_Y + 
                                  zagreb_distance +  
                                  bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 +
                                  quad1 + quad2 + quad3, data=data)


######################
#Step1: RELEVANT VARS#
######################


relevant_vars<-c("elevation", "slope", "tmp_mean_2010", "pre_mean_2010",
                 "maizemed", "river_total", "river_density",
                 "trade_route_1450_density", "trade_center_count")

#####################
#Step2: MEANS AND SD#
#####################

means<-list()
sds<-list()
for (i in relevant_vars){
  #print(i)
  item_mean<-round(mean(data[[i]], na.rm=T), 3)
  item_sd<-round(sd(data[[i]], na.rm=T), 3)
  means[[length(means) + 1]] <- item_mean 
  sds[[length(sds) + 1]] <- item_sd 
}
means<-unlist(means)
sds<-unlist(sds)

####################
#Step3: LIST MODELS#
####################


list_models<-list(lm_elevation,
                  lm_slope,
                  lm_tmp_mean_2010,
                  lm_pre_mean_2010,
                  lm_maizemed,
                  lm_river_total,
                  lm_river_density,
                  lm_trade_route_1450_density,
                  lm_trade_center_count)
                  

######################
#Step5: Column labels#
######################

column_labels= c("Elevation",
                 "Slope",
                 "\\shortstack{Avg. \\\\ Temperature}",
                 "\\shortstack{Avg. \\\\ Precipitation}",
                 "\\shortstack{Maize\\\\ Suitability}",
                 "\\shortstack{Total River\\\\ Length in km}",
                 "\\shortstack{River\\\\ Density}",
                 "\\shortstack{Trade Route\\\\ Density\\\\ 1450}",
                 "\\shortstack{No. of Trade\\\\ Centers\\\\ 1450}")



croatia_grid1=stargazer(list_models,
                       column.labels= column_labels,
                       keep = c("treat"),
                       se = list(coeftest(lm_elevation, cluster.vcov(lm_elevation, data$NAME_1))[, 2],
                                 coeftest(lm_slope, cluster.vcov(lm_slope, data$NAME_1))[, 2],
                                 coeftest(lm_tmp_mean_2010, cluster.vcov(lm_tmp_mean_2010, data$NAME_1))[, 2],
                                 coeftest(lm_pre_mean_2010, cluster.vcov(lm_pre_mean_2010, data$NAME_1))[, 2],
                                 coeftest(lm_maizemed, cluster.vcov(lm_maizemed, data$NAME_1))[, 2],
                                 coeftest(lm_river_total, cluster.vcov(lm_river_total, data$NAME_1))[, 2],
                                 coeftest(lm_river_density, cluster.vcov(lm_river_density, data$NAME_1))[, 2],
                                 coeftest(lm_trade_route_1450_density, cluster.vcov(lm_trade_route_1450_density, data$NAME_1))[, 2],
                                 coeftest(lm_trade_center_count, cluster.vcov(lm_trade_center_count, data$NAME_1))[, 2]),                       
                       covariate.labels=c("Military Territory"),
                       dep.var.caption = "Dependent variable",
                       dep.var.labels.include = F,
                       omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                       add.lines=list(c("Mean", means),
                                      c("SD", sds),
                                      c("Boundary FE", rep("Yes", 9)),
                                      c("Region FE", rep("Yes", 9))))
note_latex <- "\\multicolumn{10}{l} {\\parbox[t]{25cm}{ \\footnotesize \\textit{Notes:} 
Coefficients and municipality (opcina) clustered robust standard errors in parantheses from OLS regression. 
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.
The unit of analysis is municipality (opcina). 
See codebook in the appendix for data sources and how the variables were calculated.}}"
croatia_grid1[grepl("Note",croatia_grid1)] <- note_latex
#cat(croatia_grid1, file="croatia_grid1.tex", sep="\n")
cat(croatia_grid1, file="./Dropbox/Legacies_Project/Paper/tables/table_A3.tex", sep="\n")
